#include "coccurance.h"

// Param list contains height, width, gray scale level and
// distance between two neighbors. Default value for distance is 1 which
// stands for direct neighbor of each pixel
Coccurance::Coccurance(int h, int w, int level) :
    height(h), width(w), nGray(level)
{
    delta1 = delta2 = delta3 = 0.0;

    init();
}

/*
 *	Clear memory that stores glcm probabilities and clear all nodes
 *	in all linked lists
 */
Coccurance::~Coccurance()
{
    for(int i = 0; i < 4; i++)
    {
        for(int p = 0; p < nGray; p++)
        {
            for(int q = 0; q <= p; q++)
            {
                delete glchs[i][p][q].listnode;
            }
            delete[] glchs[i][p];
        }
        delete[] glchs[i];
        delete head[i];
    }

    delete[] glchs;
    delete[] head;
    delete[] tail;
    delete stats;
}

/*****************************
    HashNode*** glchs;
    ListNode** head;
    ListNode** tail;
    Stat* stats;

    Since symmetrical GLCM is considered, only to initialize lower half
    of 'glchs' matrix. Also one 'glchs' for each direction. Dimension of
    each 'glchs' is based on number of gray level 'nGray'. Each cell of
    'glchs' stores a listnode that contains (i, j) as well as coccurrance
    prob of(i, j)
*******************************/
void Coccurance::init()
{
	// First dimension of glchs representing 4 directions
    glchs = new HashNode** [4];
	// Array of head nodes for each direction
    head = new ListNode* [4];
	// Array of tail nodes for each direction
    tail = new ListNode* [4];
	// Pointer of statistics to store various texture features
    stats = new Stat();
	// For each direction among 0, 45, 90 and 135, index is 0,
	// 1, 2, 3 for each direction respectively. 
    for(int m = 0; m < 4; m++)
    {
		// Initialize second dimension of glchs to be number of 
		// gray levels
        glchs[m] = new HashNode* [nGray];
		// Each of the four head node points to an empty node 
        head[m] = new ListNode();
		// Let each of the tail nodes point to head node
        tail[m] = head[m];
		// Initialize 2-D hash table
        for(int i = 0; i < nGray; i++)
        {
            // Each row has increasing capacities to store coccurance
            glchs[m][i] = new HashNode[i+1];

            // Symetrical matrix is used here
            // Only lower half of matrix is used to store nodes'
            // probabilities. Initially probability of 
            for(int j = 0; j <= i; j++)
            {
                ListNode* node = new ListNode();
                node->x1 = i;
                node->x2 = j;
                node->next = 0;
                glchs[m][i][j].prob = 0.0;
                glchs[m][i][j].listnode = node;
            }
        }
		// After initialization, have the next node of each head node
		// point to the first node in glchs 
        head[m]->next = glchs[m][0][0].listnode;
    }   
}

/*
 *  Reset all variables, so each cell within 'glchs' will not
 *  point to any other nodes, and set all probabilities value
*	to be 0. Then reset head and tail of all
 *  four directions.
 */
void Coccurance::reset()
{
    // For each direction
    for(int m = 0; m < 4; m++)
    {
        // For each value of number of gray levels
        for(int i = 0; i < nGray; i++)
        {
            // Symmetrical matrix is used here. Set value of probability
			// and the next listnode to zero(null)
            for(int j = 0; j <= i; j++)
            {
                glchs[m][i][j].listnode->next = 0;
                glchs[m][i][j].prob = 0.0;
            }
        }
        // Points head's next node to null
        head[m]->next = 0;
        // Points tail node to head
        tail[m] = head[m];
    }
    // Reset each feature value to 0's
    stats->con = 0.0;
    stats->dis = 0.0;
    stats->uni = 0.0;
    stats->idm = 0.0;
    stats->ent = 0.0;
    stats->cor = 0.0;
    stats->inv = 0.0;
    stats->var = 0.0;
}


/*  
 *	Passing address of first element pointed by 'data'
 *  as well as distance 'd' between 2 neighbor pixels.
 *  Then build link list for each direction and calcuate
 *  texture features.
 */
void Coccurance::setParam(int *data, int d)
{
    // Prerequisites for calculation of coccurrance
    if(data == 0 || d >= height || d >= width || d <= 0) return;

    // Total number of coccuring pairs
    // Symmetrical matrics is considered
    // delta1 is for 0
    // delta2 is for 45 and 135
    // delta3 is for 90
    delta1 = 1.0 / ((float) height*(width-d));
    delta2 = 1.0 / ((float) (width-d)*(height-d));
    delta3 = 1.0 / ((float) (height-d)*width);

    // build coccurance matrix for int *d
    build_list(data, d);
    // calculate texture features
    calc_stat();
}

/*  
 *	Build coccurrance for two neighbor pixels at distance 'd'
 *  There are four directons to be considered neighbor:
 *  0, 45, 90 and 135 degree. 'data' is enclosed in region of
 *  size (height x width). Calculation are separated into four
 *  regions. Within each region, angle/directions are different
 *  with regard to co-occurrence.
 */
void Coccurance::build_list(int *data, int d)
{
    // i - height
    // j - width
    int i, j;
    // For each row from d to height-1
    for(i = d; i < height; i++)
    {
        // For each column from d to width-d
        for(j = d; j < width-d; j++)
        {
            // Within this region, all four direction's neighbor
            // are available for co-occurrence calculation
            find_insert(*(data+i*width+j), *(data+i*width+j+d), 0);
            find_insert(*(data+i*width+j), *(data+(i-d)*width+j+d), 45);
            find_insert(*(data+i*width+j), *(data+(i-d)*width+j), 90);
            find_insert(*(data+i*width+j), *(data+(i-d)*width+j-d), 135);
        }
    }
    // For each row from 0 to d-1
    for(j = 0; j < d; j++)
    {
        // For each column from d to height-1
        for(i = d; i < height; i++)
        {
            // Within this region, only 0, 45, and 90 degree angle
            // are available for co-occurrence calculation
            find_insert(*(data+i*width+j), *(data+i*width+j+d), 0);
            find_insert(*(data+i*width+j), *(data+(i-d)*width+j+d), 45);
            find_insert(*(data+i*width+j), *(data+(i-d)*width+j), 90);
        }
    }
    // For each row from 0 to d-1
    for(i = 0; i < d; i++)
    {
        // For each column from 0 to width-d-1
        for(j = 0; j < width-d; j++)
        {
            // Within this region only 0 degree angle is available
            // for co-occurrence calculation
            find_insert(*(data+i*width+j), *(data+i*width+j+d), 0);
        }
    }
    // For each row from width-d to width-1
    for(j = width-d; j < width; j++)
    {
        // For each column from d to height-1
        for(i = d; i < height; i++)
        {
			// Within this region, only 90 and 135 degree co-occurrence 
			// are available for calculation.
            find_insert(*(data+i*width+j), *(data+(i-d)*width+j), 90);
            find_insert(*(data+i*width+j), *(data+(i-d)*width+j-d), 135);
        }
    }
}

/*
 *  Given a pair (x1, x2) at 'angle', find out the probabilities
 *  between them, known as cooocurance. If it is 0, then create a new
 *  node and append to the end of linked list. If it is not 0,
 *  simply update the probabilities of existed node in glchs matrix.
 *
 */
void Coccurance::find_insert(int x1, int x2, int angle)
{
    // Convert angle to index. angle values are 0, 45, 90, 135
    // index values are 0, 1, 2, 3.
    int index = angle / 45;

    // Quantize value down to 'ng' level
    x1 = (int) x1*nGray/256;
    x2 = (int) x2*nGray/256;

    // Force x1 >= x2 for lower half matrix
    if(x1 < x2)
    {
        int temp = x1;
        x1 = x2;
        x2 = temp;
    }
    // Get corresponding co-occurrence probabilities as well
    // as list node from linked list based on index, x1 and x2
    float prob = glchs[index][x1][x2].prob;
    ListNode* n = glchs[index][x1][x2].listnode;

    // If prob==0.0, then append it to the end of sequence
    // pointed by tail pointer. This is because that this
    // pair of pixels/elements has not been encountered
    if(prob == 0.0)
    {
        tail[index]->next = n;
        tail[index] = n;
    }
    // Depending on different angle, update probability of entry
    // in glchs by adding different delta value.
    if(angle == 0)
    {
        glchs[index][x1][x2].prob = prob + delta1;
    }
    else if(angle == 45)
    {
        glchs[index][x1][x2].prob = prob + delta2;
    }
    else if(angle == 90)
    {
        glchs[index][x1][x2].prob = prob + delta3;
    }
    else if(angle == 135)
    {
        glchs[index][x1][x2].prob = prob + delta2;
    }
}

/*
 *	Calculate statistics for current WSIZE x WSIZE
 *	block. Starting from the head node. Final value
 *	will be average of all four directions.
 */
void Coccurance::calc_stat()
{
	// Current node in the linked list
    ListNode* t;

    // Contrast, Uniformity, Dissimilarity, Entropy
    // Inverse Difference Moment, Correlation,
    // Inverse Difference and covariance
    float con=0.0, uni=0.0, dis=0.0, idm=0.0;
    float ent=0.0, cor=0.0, inv=0.0, var=0.0;
    float avg=0.0, variance=0.0, prob=0.0;
	// Location of index in glchs
    int x1 = 0, x2 = 0;
    // Go through each angle, 0, 45, 90 and 135
    for(int i = 0; i < 4; i++)
    {
        // Obtain first listnode of each direction
        t = head[i]->next;
        // Reset average and covariance values
        avg = variance = 0.0;
        // Get average and covariance for current linked list
		// Then add to final variance result which is one of 
		// the texture features
        if(t != 0)
        {
            getAvgVar(i, t, avg, variance);
            var += variance;
        }
        // Iterate through the linked list of current angle/direction
        while(t != 0)
        {
			// Get current location index x1 and x2 pointed by 't' node
            x1 = t->x1;
            x2 = t->x2;
			// Get current probability pointed by 't' node
            prob = glchs[i][x1][x2].prob;

			///////////////////////////////////////////////////////
			///////////////////////////////////////////////////////
			// Harralick texture feature calculation 
            con += prob * (x1-x2)*(x1-x2);
            dis += prob * abs(float (x1-x2));
            uni += prob * prob;
            idm += prob / (1+(x1-x2)*(x1-x2));
            ent += prob * (log(prob));
            cor += (x1-avg)*(x2-avg) * prob/variance;
            inv += prob / (1+x1-x2);
			///////////////////////////////////////////////////////
			///////////////////////////////////////////////////////

			// Move to the next node in the link list
            t = t->next;
        }
    }

    // Obtain average across all four directions
    stats->con = con / 4.0;
    stats->uni = uni / 4.0;
    stats->dis = dis / 4.0;
    stats->idm = idm / 4.0;
    stats->ent = -1*ent / 4.0;
    stats->cor = cor / 4.0;
    stats->inv = inv / 4.0;
    stats->var = var / 4.0;
 }


/*
 *	Return calculated feature value based on 'func' which is 0-based
 */
float Coccurance::avg_stat(int func)
{
	// Value of retrieved feature value
    float value = 0.0;
	// Based on input 'func' value which is among [0..7]
    switch (func)
    {
    case 0 :    // Contrast
        value = stats->con;
        break;
    case 1 :    // Dissimilarity
        value = stats->dis;
        break;
    case 2 :    // Uniformity
        value = stats->uni;
        break;
    case 3 :    // Inverse difference moment
        value = stats->idm;
        break;
    case 4 :    // Entropy
        value = stats->ent;
        break;
    case 5 :    // Correlation
        value = stats->cor;
        break;
    case 6 :    // Inverse difference
        value = stats->inv;
        break;
    case 7 :    // Covariance
        value = stats->var;
        break;
    default:
        break;
    }

    return value;
}

/*
 *   Get average and covariance for linked list. The node
 *   is the head node of each direction indicated by index
 */
void Coccurance::getAvgVar(int index, ListNode *node, float &avg, float &var)
{
    int x1 = 0,  x2 = 0;
    float prob = 0.0;
	// Node is the head node each linked list
    ListNode *t = node;

    // Calculate average of the linked list
    while(t != 0)
    {
		// Get current node's location index x1 and x2, plus 
		// probability value
        x1 = t->x1;
        x2 = t->x2;
		// prob has been quantized
        prob = glchs[index][x1][x2].prob;
		// Add to average value 
        avg += x1 * prob;
		// Move to next node
        t = t->next;
    }
    // Reset t to  point to start of linked list
    t = node;
    // Calculate covariance of the linked list
    while(t != 0)
    {
        x1 = t->x1;
        x2 = t->x2;
        prob = glchs[index][x1][x2].prob;
		// variance calculation
        var += (x1 - avg) * (x1 - avg) * prob;
		// Move to next node
        t = t->next;
    }
	// Reset t to null
    t = 0;
}
